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First-principles density functional theory (DFT) oxygen chemical shift tensors were calculated for 
A(B,B')03 perovskite alloys Pb(Zri/2Tii/2)03 (PZT) and Pb(Mgi/3Nb2/3)03 (PMN). Quantum 
chemistry methods for embedded clusters and the GIPAW method [C. J. Pickard and F. Mauri, Phys. 
Rev. B 63 245101 (2001)] for periodic boundary conditions were used. Results from both methods 
are in good agreement for PZT and prototypical perovskites. PMN results were obtained using only 
GIPAW. Both isotropic Jiso and axial (5ax chemical shifts were found to vary approximately linearly 
as a function of the nearest-distance transition-metal/oxygen bond length, r^. Using these results, 
we argue against Ti clustering in PZT, as conjectured from recent ^^O NMR magic-angle-spinning 
measurements. Our findings indicate that ^^O NMR measurements, coupled with first-principles 
calculations, can be an important probe of local structure in complex perovskite solid solutions. 

PACS numbers: 76.60.Cq 77.84.-s 
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I. INTRODUCTION 

High performance solid solution ferroelectrics, based 
on the ideal ABO3 perovskite structure, are widely used 
in technological applications such as ultrasonic trans- 
ducers, sensors, actuators, and thin film applications^^ 
The strong electromechanical coupling in these materials 
is related to a balance of competing instabilities, such 
as cation off-centerings and oxygen octahedral rotations. 
The perovskite structure offers many ways to fine tune 
these interactions through chemical substitutions, such 
as alloying on the A and/or B sites, and through epitax- 
ial control in layered and thin film geometries.— 

Solid state nuclear magnetic resonance (NMR) has in- 
creasingly been used to study the local structure and 
dynamics of these complex perovskites i^"— NMR spectra 
of a target nucleus are largely determined by the cou- 
pling of its magnetic dipole and electric quadrupolc mo- 
ments with the local magnetic field and electric field gra- 
dient, respectively. The interpretation of NMR spectra 
in complex solid solution perovskites is complicated by 
the presence of broad spectral features due to disorder. 
First-principles calculations of electric field gradient^ and 
chemical shielding^ tensors can play an important role 
interpreting NMR spectra in these materials. 

Previously we demonstrated a near linear dependence 
of the chemical shielding tensor a on the oxygen near- 
est neighbor B-0 bond distance Vg for the prototypical 
perovskites BaTiOg (BT), SrTiOg (ST), PbTiOg (PT), 
and PbZrOs (PZ)4 The linear dependence was shown to 
arise from large paramagnetic contributions to a^ and ay 
principal values (our convention identifies the "z" princi- 
pal axis as that most nearly parallel to the B-O-B bond 
direction), due to virtual transitions between 0(2p) and 
unoccupied B(nd) states. This linear variation is con- 
firmed here for two complex perovskite solid solutions, 
PZT and PMN. 



First-principles chemical shielding calculations have 
traditionally been done with the embedded cluster ap- 
proach, using standard quantum chemistry methods^i"— 
with gaussian type orbitals (GTO). More recently, the 
planewave based GIPAW method with PBC has provided 
an alternative approach. Relatively few calculations for 
transition metal oxides have been reported using either 
technique. Here we use complementary calculations with 
both methods to cross validate convergence with respect 
to cluster size and termination effects, basis sets, and the 
accuracy of pseudopotentials (PSPs). 

Quantum-chemistry methods can calculate chemical 
shielding tensors for embedded clusters, using a range 
of approximations, from Hartree Fock and density func- 
tional theory (DFT) with semilocal or hybrid exchange- 
correlation functionals, to explicitly correlated meth- 
ods such Moller-Plesset perturbation theory and coupled 
cluster approaches.—"— The principal difficulties with 
the embedded cluster approach are controlling size and 
basis set convergence. Size effects can be monitored by 
studying increasingly larger clusters. Long-range elec- 
trostatic interactions can be handled by embedding the 
cluster in large arrays of point charges and eliminating 
depolarizing fields, as described in Ref. [^ Achieving the 
basis set limit can be problematic in some cases, because 
atom-centered GTO's do not form a complete orthogo- 
nal basis. Nevertheless, basis set convergence is generally 
well controlled through the use of standardized GTO ba- 
sis sets.— The cluster approach becomes inefficient, how- 
ever, for complex systems, because separate cluster cal- 
culations are usually required for each inequivalent target 
atom. 

GIPAW calculations are naturally adapted to or- 
dered crystalline solids, since size convergence is effec- 
tively achieved by using primitive unit cells with well- 
converged /c-point quadrature grids for Brillouin zone in- 
tegrations. Disordered solids can be treated using super- 
cells. Planewaves form a complete basis, so convergence 



to the basis set limit is straightforward. The method ap- 
pHes the projector-augmented-wave (PAW) reconstruc- 
tion^ to a conventional PSP calculation to obtain all- 
electron valence wave functions, which are required for 
accurate calculations of the chemical shielding. There 
are two principle issues with the GIPAW method. GI- 
PAW PSPs are more difficult to construct than standard 
norm-conserving PSPs. To achieve good transferability, 
they may require multiple partial wave channels and large 
planewave cutoff energies for some target atoms. The 
construction of the PAW atomic-like augmentation basis 
also requires care. To date, GIPAW calculations have 
been carried out for only a limited number of transition 
metal oxide systemsji^i^^ 

PZT is a homovalent mixture of Ti'*^ and Zr"*^ transi- 
tion metal cations, while PMN is a heterovalent 2:1 mix- 
ture of the Nb^"^ transition metal cation and the Mg^+ 
alkaline earth cation. Both embedded cluster and GI- 
PAW calculations were carried out for PZT, while PMN 
results were obtained only with the GIPAW method, us- 
ing PBC. First-principles relaxed structural models were 
used to simulate PZT and PMN structures. 

The remainder of the paper is organized as follows. 
The theoretical approaches are described in Sec. HH Re- 
sults and discussion are presented in Sees. Hill and IIV| 
respectively. We summarize and conclude in Sec. [V] 



II. THEORETICAL METHODS 

The chemical shielding tensor a determines the total 
magnetic field at an atomic nucleus. 



B = (l-a)Be 



(1) 



where Boxt is the external field. For the systems consid- 
ered here, a is calculated using embedded cluster and 
GIPAW-PBC methods. The symmetrioi^i^o a tensor 
is determined by its principal axis components, with 
isotropic and anisotropic parts, conventionally defined 



CTax = g (20-2 - 



~ ^y) - 2 



= iTra 

2 (O'z ^ 0-iso) 



(2) 



As mentioned, our convention for the perovskite struc- 
ture identifies the "z" principal axis as that most nearly 
parallel to the B-O-B bond direction. NMR measure- 
ments of a are usually reported with respect to a refer- 
ence material, where the chemical shift tensor 5 is defined 



This yields the relation 






(4) 



Rather than regarding the slope m as an independent 
fitting parameter, we set tti = 1, which leaves cr'^j as 
the only independent parameter. Allowing m to vary 
yields fits of similar overall quality, but somewhat dis- 
torts chemical shift differences between inequivalent oxy- 
gen sites within the same material. Constraining m = 1 
allows better cancellation of errors, due to systematic ef- 
fects such as the choice of methodological approach or 
DFT exchange-correlation functional. 



A. Embedded cluster calculations 

A detailed discussion of this approach is given in 
Ref. 1^. We briefiy summarize some of the key fea- 
tures of this method. With central O atoms, embedded 
clusters consisted of either 21 "quantum" (QM) atoms, 
(A4B2O15)"", or 65 QM atoms, (A4Bio05i)^^". In 
these clusters, all cation atoms are fully coordinated with 
QM O atoms. The QM clusters arc further surrounded 
by a large array of point charges, which reproduce the 
Madelung potential in the QM region»^ To alleviate arti- 
ficial polarization of boundary 0(2p) states, the nearest- 
neighbor (nn) and the next-nearest-neighbor cation point 
charges of boundary O atoms are replaced by "empty" 
PSPs (ePSPs)<^i2i Finally, in non-ccntrosymmetric clus- 
ters, depolarizing electric fields are removed by applying 
an external electric field<^ 

Calculations were performed with the GAUSSIAN com- 
putational package ^iii and the chemical shielding ten- 
sor was determined using the continuous set of gauge 
transformations (CSGT) mcthodi^i^ Calculations were 
done using the DFT hybrid B3LYP— , as well as general- 
ized gradient approximations (GGA), using the PW91— 
and PBE^ forms. Douglas-Kroll-Hess 2nd order scalar- 
rclativistic calculations were performed on selected sys- 
tems. Atom-centered GTO basis functions were asso- 
ciated with all the QM atoms. All-electron treatments 
were used for the O and Ti atoms, while the other QM 
atoms were represented using scalar-relativistic small 
core (scalar-RSC) PSPs [also called effective core poten- 
tials (ECPs)]. The well- converged GTO basis sets and 
ECPs used for these calculations are described in Ref. [^ 
and were taken from the EMSL websiteJ^ 



B. GIPAW calculations 



6 = -{a - CTiof), 



(3) 



with corresponding definitions to those in Eq. [51 The 
experimental reference for oxygen is liquid H2O. 

The theoretical oxygen reference value cr'^f is deter- 
mined from a linear regression of S^^^ ' versus (y-,s^™^ J^i^ 



Calculations with PBC used the GIPAW functional- 
ity of the QuantumESPRESSO (QE) codei^ These calcu- 
lations are done in two steps. A standard ground state 
norm-conserving PSP calculation is first performed. This 
is followed by a linear response calculation in the pres- 
ence of an external magnetic field. The linear response 



TABLE I: Pseudopotential construction parameters, see text for description. 





reference state 


rc(au) 


reference energies (Ry) 


core correction radius (au) 





2s^2p^,5rf,5s,5p 


1.2,1.2,1.5 


*,*,0. 10,0. 10,0. 10 


0.34 


Ti 


3s^,3p^, 3d°,4s,4p,4d 


0.9,0.9,0.9 


* 


* 


*,-2. 73,0.50,0.10 


0.53 


Pb 


Qs^fip°M^°,7s,7p,6d 


2.0,2.4,1.0 


* 


* 


*,-1.00,-0.20,-1.80 


0.89 


Zr 


4s^4p6,4d°, 5s, 5p,5d 


1.0,1.2,1.4 


* 


* 


*,-2.45,-2.00,-1.37 


0.80 


Sr 


4s2,4p6,4d°, 5s, 5p 


1.2,1.3,1.7 


* 


* 


-1.00,-0.72 


0.88 


Ba 


bs"^, 5p^, 5 d, 6s, 6p 


1.5,1.7,2.0 


* 


* 


-0.95,-0.90,-1.50 


1.19 


Mg 


2s'^,2p^,3d,3s,3p 


0.6,0.6,1.5 


* 


* 


0.3,-1.3,0.1 


0.30 


Nb 


4s2,4p6,4d°, 5s, 5p,5rf 


1.0,1.1,1.3 


* 


* 


*,-3.30,-2. 70,-1.27 


0.80 


K 


3s2,3p6,3d°,^s,^p 


1.4,1.5,1.5 


* 


* 


*,-1.50,-0.80 


0.50 



calculation uses all-electron like valence wave functions, 
which are represented by planewaves, modified near the 
nuclei by atomic-like PAW augmentation basis functions. 
The PAW basis functions are constructed when the PSP 
is generated, as further discussed below. Calculations 
reported below used well converged Monkhorst-Pack^i 
Brillouin zone fc-point sampling, e.g. 6x6x6 for the pro- 
totypical perovskites ST, BT, and PT. Unless otherwise 
specified, the PBE GGA functional was used in all QE 
calculations. 



All PSPs were constructed using the "Idl" PSP gener- 
ation code (distributed with QE). All norm-conserving 
PSPs were scalar-relativistic Troullier-Martins^ type. 
Table |T] shows the construction parameters used for all 
PSPs employed in this paper. The local channel is in- 
dicated by boldface type. States that are in italics were 
generated using a Hamann type^ reference state. The 
re's correspond to the s, p, and d channels respectively. 
The next column shows the PSP reference energies for 
the s, p, and d channels. The symbol '*' indicates that 
the corresponding all-electron eigenvalue was used for 
this state. It should be noted that these parameters 
are for the PSP projectors used in the ground state self- 
consistent total-energy calculations. Parallel to the PSP 
construction, a second atomic calculation is performed to 
generate the required GIPAW augmentation basis func- 
tions, consisting of all-electron and pseudo partial wave 
radial functions. The number of GIPAW angular mo- 
mentum channels was the same as for the PSPs. For 
all except Pb, O and K, the augmentation basis func- 
tions were generated using the corresponding all-electron 
eigenvalues (n.6., the Ti^"*" Ad state is bound, for exam- 
ple), rather than the values in Table HI For Pb, O and 
K, the values in Table |T] were used, except for the Pb 7s, 
which used the all-electron eigenvalue. The final column 
shows the non-linear core-correction radius for each po- 
tential. A conservative 350 Ry energy cutoff was used. 
This could have been reduced by using larger re's for the 
metal PSPs. The 350 Ry energy cutoff, while high, was 
easily tractable for all systems studied in this paper. This 
cutoff yields chemical shieldings to within about 1 ppm, 
as indicated by test calculations with other settings. 



TABLE IL Comparison of chemical shielding results, u sing 
embedded clusters and GIPAW-PBC from QE and Ref. [3 
The principle values of the oxygen chemical shielding tensor 
are presented for the Ti02 molecule and rutile, cubic ST, cu- 
bic BT, and tetragonal PT. All values are from GGA calcula- 
tions. Embedded cluster results are labeled as C-n, where n is 
is the number of QM atoms in the cluster. The GIPAW-PBC 
and C-65 calculations were done with relativistic PBE, while 
the C-21 calculations were done with non-relativistic PW9L 



CTx 


(Jy 


Oz 


^iso 


TiOa (molecule) 








Gaussian -1803 


-801 


148 


-819 


GIPAW -1826 


-811 


146 


-830 


Ti02 (rutile) 








C-99 -499 


-380 


-306 


-395 


GIPAW -483 


-380 


-296 


-386 


SrTiOg (cubic) 








C-21 -353 


-353 


46 


-220 


C-65 -403 


-403 


27 


-260 


GIPAW -429 


-429 


7 


-284 / -287" 


BaTiOs (cubic) 








C-21 -414 


-414 


49 


-260 


C-65 -483 


-483 


22 


-315 


GIPAW -529 


-529 


-31 


-363 / -379* 


PbTiOs-axial (P4mm) 






C-21 -562 


-562 


123 


-334 


C-65 -599 


-599 


67 


-377 


GIPAW -630 


-630 


77 


-394 


PbTiOs-equatorial (P4mm) 






C-21 -286 


-228 


-32 


-182 


C-65 -365 


-277 


-35 


-226 


GIPAW -398 


-284 


-23 


-235 


"Ref. 


li determined Uiso 


using the 


experimental 


structure 


''Ref. 


la determined Ui^o 


using the relaxed structure 



C. Comparison of methods 

Table HIl shows the comparison between the GIPAW 
and cluster approach for the Ti02 molecule, the rutile 
crystalline solid, as well as some results for prototypical 



TABLE III: The derived theoretical oxygen reference value, cr'cfi is used to determine theoretical isotropic chemical shifts 
from the corresponding calculated isotropic chemical shielding values. Experimental chemical shifts are shown for comparison. 
The values of aj^f a-re shown in the third row from the bottom of this table. The rms error and maximum deviation of the 
calculated shifts, compared to experiment, are also shown. Calculated cluster results (21 and 65 QM atoms) with B3LYP and 
PW91 exchange-correlation are shown together with GIPAW-PBC with PBE exchange-correlation. (PZ experimental chemical 
shift site assignments were corrected in Ref. [3 and are used here.) 





Expt" 


Cluster 

B3LYP PW91 

21 65 21 65 


present ^ 


PBC-GIPAW 

present 


other 


ST (cubic) 


467 ±5 


491 


477 


494 


491 


491 


496 


503 '' 


BT (cubic) 


546 ±5 


536 


537 


534 


547 


570 


575 


595 " 



BT (P4nim) 






















Oax 


570 ±5 


591 




571 






579 


573 '^ / 611 ^ 




Oeq 


520 ±5 


515 




516 






535 


563 "^ / 531 " 




A(Oax - Oeq) 


50 


76 




56 






44 


10/80 


PT (P4mm) 






















Oax 


647 ±2 


644 


644 


608 


609 


601 


606 






Oeq 


443 ±2 


449 


445 


457 


455 


442 


447 






A(Oax - Oeq) 


204 


195 


199 


152 


154 


159 


159 




PZ (Pbani) 






















01-4g 


365 ±2 


357 




367 






355 






01'-4g 


351 ±2 


346 




356 






336 






02-81 


356 ±2 


355 




364 






349 






03-4f 


329 ±2 


324 




336 






309 






04-4e 


408 ±2 


392 




399 






415 






A(04-01) 


43 


35 




31 






60 






A(04 - 01') 


57 


47 




43 






79 






A04 - 02) 


52 


38 




35 






66 






A(04 - 03) 


79 


68 




63 






105 





„th. 


288 " 


293 


238 


275 


217 


207 


212 


216 


RMS error 




12 


7 


16 


23 


28 


20 


26 


Max Dev ^ 




24 


10 


39 


38 


46 


41 





"ST and BT experimental chemical shift values are from Ref. [y and the PT and PZ experimental chemical shift values are from Ref. |j. 

'Ref. [13 determined Uiso using the experimental structure 

'^Ref. [13 determined ai^io using the relaxed structure 

"^The ^'^O experimental chemical shielding reference is liquid water, o""j^*°'' = 287.5 ppm24 

«Max Dev = max(abs(5'h<=o--y _ ^cxpt.-,) 

^Used (T*^j determined from linear regression on the same systems as the 65 QM atom clustcr-PW91 results (sec text). 



perovskites. Where available, published GIPAW results 
from Ref. [l8| are also given for comparison. 

Ti02 molecule calculations were for a Ti-0 bond length 
of 1.651 A and a 0-Ti-O angle of 114.2°. With PBC- 
GIPAW, a 22x24x28 Bohr superccU was used. The 
Gaussian calculation was for the isolated molecule and 
used a basis set of cc-pwCVQZ-DK and IGLO-III for 
Ti and O, respectively. The experimental rutile struc- 
ture was used for sohd TiOa-^ The PBC-GIPAW cal- 
culation used a fc -point sampling of 4x4x6. The cluster 



method used a 077Ti22 QM cluster with cc-pwCVTZ-DK 
and 6-311G(2d,p) basis sets for the inner 3 and outer 19 
Ti atoms, respectively; IGLO-III and IGLO-II basis sets 
were used for the inner 15 and outer 52 O atoms, respec- 
tively. (Embedding with ePSPs and point charges was 
done as described above.) 

Calculations for the prototypical perovskites in Table 
iniused the experimental structures as described in ReL|^. 
PBC-GIPAW used a fc-point sampling of 6x6x6. Em- 
bedded clusters were constructed as above (Section lll Ap . 



and results are given for 21 and 65 QM atom clusters. 
The 21 QM atom cluster results are from non-relativistic 
calculations using the PW91 DFT functional. The 65 
QM atom clusters are from scalar relativistic PBE calcu- 
lations. Differences between PW91 and PBE (not shown 
in the Table) and relativistic effects are small, as ex- 
pected. Tests on the 65 QM atom ST cluster show that 
non-relativistic PBE produced cr^.y and a^ values that 
are ci 6 and 1.5 ppm more positive, respectively, than 
PW91. Adding scalar-relativistic effects changed a^^y 
and Oz by ~ -1-12 and -1-3 ppm, respectively, indepen- 
dent of the GGA functional. These effects are negligible 
for the corresponding chemical shifts, as expected, due 
to cancellation of errors. 

The ST and PT 65 QM atom relativistic PBE cr^,j^ and 
Oz are in good agreement with PBC values, differing at 
most by ~ 33 and 20 ppm respectively. Isotropic values 
for both systems arc within 24 ppm of PBC-GIPAW. 
A larger discrepancy is seen in BT, where the cluster 
values are rigidly shifted by about -1-50 ppm. Size effects 
between the 21 and 65 QM atom results are evident in 
the table. The 21 QM atom values are more shielded [i.e., 
more positive (see Eq. ([T|))] than either the 65 QM atom 
values or the PBC methods. Rigidly shifting all 21 QM 
atom brings these into better agreement, indicating good 
cancellation of errors. This is evident in the chemical 
shifts shown in Table IIIII in the next section, where 21 
and 65 QM atom clusters are seen to give nearly identical 
chemicals shifts. 

These results demonstrate that the embedded cluster 
and PBC-GIPAW approaches produce comparable agree- 
ment with measured isotropic chemical shifts, regardless 
of cluster size and methodology. Differences due to clus- 
ter size, DFT functionals, relativistic effects, and PSPs 
largely cancel in the chemical shifts, i.e., they are ab- 
sorbed in the constant chemical shielding reference value 
th. 

"ref- 



ill. RESULTS 

In this section, we first describe the calculation of 
the theoretical oxygen chemical shielding reference af^^ 
[Eq. (U))]. As mentioned, this is done using linear regres- 
sion of the calculated chemical shieldings with the corre- 
sponding measured chemical shifts, where available.— i^ 
We next present calculated ^^O chemical shifts for two 
perovskite-based B-site alloys, PZT and PMN, using the 
derived values of a^^^ . 



A. Determination of the theoretical ^^O chemical 
shielding reference 

A linear regression was separately evaluated for the 
PBC and embedded cluster calculations J^i^^ Addition- 
ally for the clusters, separate regressions were performed 
for different cluster sizes and DFT functionals. Table Hill 



shows derived af^^ values for each case, along with the 
rms error and maximum deviation in the isotropic chem- 
ical shifts. The table compares the resulting calculated 
isotropic chemical shifts to the measured values. As in 
Ref. |8J, experimental structures were used for all sys- 
tems except PZ. For PZ, experimental lattice parame- 
ters from neutron scattering measurements were used to- 
gether with internal coordinates determined from first- 
principles calculations)^ 

Calculated cluster results (21 an 65 QM atoms) 
with B3LYP and PW91 exchange-correlation are 
shown together with GIPAW-PBC with PBE exchange- 
correlation. The last column shows other GIPAW results, 
where available. For the 65 QM clusters, crJJ^f is derived 
from a more limited set of calculations, as shown in the 
Table. For comparison, the effect of this to be only a 
few ppm. B3LYP results are seen to give slightly better 
agreement with experiment. Using the same exchange- 
correlation treatment, both the 21 and 65 atom QM clus- 
ters are seen to give nearly identical chemicals shifts. The 
values of aJJ^j in Table IIIII are used below to determine 
the theoretical chemical shifts for PZT and PMN in Ta- 
bles IIVI and |Vl Given the small difference between the 
PW91 and PBE GGA functionals and cancellation effects 
in chemical shifts, both PW91 and PBE arc labeled as 
GGA in all further results below. 



B. Results for V\,{T.vx-^'YU)Oo^ (PZT) and 

Pb(Mgi/3Nb2/3)03 (PMN) 

PZT is a homovalent mixture of Ti^+ and Zr^"*" tran- 
sition metal cations, while Pb(Mgi/3Nb2/3)03 (PMN) is 
a heterovalcnt 2:1 mixture of the Nb^"*" transition metal 
cation and the Mg^+ alkaline earth cation. Both em- 
bedded cluster and GIPAW calculations were carried out 
for PZT, while PMN results were obtained only with the 
GIPAW method, using PBC. 

Disordered PZT 50/50 was modeled using first- 
principles relaxed structures^, with different imposed 
B-site ordering and symmetry: i) [001] ordering par- 
allel to the ferroelectric polarization, with tetrag- 
onal P4mm symmetry (a x a x 2c; c/a = 1.045); ii) 
[001] ordering perpendicular to the [100] ferroelec- 
tric polarization, with orthorhombic P2mm symmetry 
(a' X a X 2a; a'/a = 1.04); iii) [111] ordering (rocksalt B 
sublattice) parallel to the ferroelectric polarization, with 
R3m symmetry. Theoretical PZT isotropic, axial, and 
anisotropic chemical shifts are summarized in Table IIVI 
[As mentioned, the axial and anisotropic components in 
Eq. ([2]) were calculated using the convention that the 
"z" principal axis is that most nearly parallel to the 
B-O-B bond direction of the target O atom.] For com- 
parison, results for ST, BT, PT and PZ are reproduced 
from Table IIIII For each inequivalent target O atom 
in the above structural models, the corresponding 21 
QM atom embedded cluster was used to calculate the 
chemical shielding tensor with both B3LYP and GGA 



TABLE IV: Calculated oxygen isotropic, axial and anisotropic components (ppm) of the chemical shift tensor for three PZT 
50/50 structural models from Ref. 0. The notation B~0~B indicates O atoms with two equidistant nn B atoms, and B-0~B 
indicates an O atom with one short and one long nn B bond. For cluster-GGA results, numbers in parenthesis show the 
difference with GIPAW. For cluster-B3LYP results, numbers in square brackets show the difference with cluster-GGA values. 
For cases where 5aniso = by symmetry, this is indicated by a dash. 







5iso 






(^ax 






^aniso 






Cluster 


GIPAW 


Cluster 


GIPAW 


Cluster 


GIPAW 




B3LYP 


GGA 


GGA 


B3LYP 


GGA 


GGA 


B3LYP 


GGA 


GGA 


ST (cubic) 


491 [-3] 


494 (-2) 


496 


-144 [-11] 


-133 (12) 


-145 


_ 


_ 


_ 


BT (cubic) 


536 [2] 


534 (-41) 


575 


-170 [-16] 


-154 (12) 


-166 


_ 


_ 


_ 


BT (P4mni) 




















Oax (Ti-O-Ti) 


591 [20] 


571 (-8) 


579 


-207 [-24] 


-183 (3) 


-186 


- 


" 


- 


Oeq (Ti-O-Ti) 


515 [-1] 


516 (-19) 


535 


-155 [-14] 


-140 (7) 


-147 


-24 [2] 


-26 (15) 


-41 


PT (P4nim) 




















Oax (Ti-O-Ti) 


644 [35] 


608 (2) 


606 


-257 [-29] 


-228 (8) 


-236 


- 


- 


- 


Ocq (Ti-O-Ti) 


449 -8] 


457 (10) 


447 


-85 [-10] 


-75 (31) 


-106 


-22 7] 


-29 (28) 


-57 


PZ(Pbam) 




















01-4g (Z-O-Z) 


357 [-10] 


367 (12) 


355 


-75 [-5] 


-70 (20) 


-90 


-33 [1] 


-34 (12) 


-46 


Or-4g (Z-O-Z) 


346 [-10] 


356 (20) 


336 


-65 [-4] 


-61 (21) 


-81 


-14 [2] 


-16 (11) 


-27 


02-8i (Z-O-Z) 


355 [-9] 


364 (16) 


349 


-81 [-4] 


-77 (14) 


-91 


-4[1] 


-5(2) 


-7 


03-4f (Z-O-Z) 


324 [-11] 


336 (26) 


309 


-47 [-4] 


-43 (9) 


-52 


-15 [-1] 


-14 (3) 


-17 


04-4e (Z-O-Z) 


392 [-6] 


399 (-16) 


415 


-137 [-7] 


-130 (36) 


-166 


-14 [1] 


-15 (11) 


-26 


PZT (P4mni) 




















Ol(Zr-O-Zr) 


352 [-9] 


361 (33) 


328 


-72 [-6] 


-66 (20) 


-85 


-16 [-4] 


-12 (-4) 


-8 


02(Zr-0-Ti) 


398 [7] 


391" (-19) 


410 


-126 [-8] 


-118 (17) 


-135 




- 


- 


03(Ti-0-Ti) 


418 [-18] 


436 (9) 


427 


-76 [-11] 


-65 (21) 


-86 


-59 [4] 


-63 (38) 


-101 


04 (Ti-O-Zr) 


679 [40] 


639 (13) 


626 


-266 [-34] 


-232 (14) 


-246 




- 


- 


PZT (P2mni) 




















Ol(Zr-O-Zr) 


358 [-9] 


367 (19) 


348 


-80 [-6] 


-74 (22) 


-96 


-36 [2] 


-38 (20) 


-59 


02(Zr-0-Zr) 


401 [-4] 


405 (-1) 


407 


-124 [-6] 


-118 (18) 


-137 


-34 [2] 


-36 (9) 


-44 


03(Ti-0-Ti) 


408 [-11] 


419 (24) 


395 


-65 [-8] 


-57 (17) 


-74 


-6 [5] 


-1 (26) 


-27 


04(Ti-0-Ti) 


668 [39] 


629 (12) 


617 


-264 [-32] 


-232 (6) 


-238 


-49 [3] 


-46 (8) 


-54 


05(Ti-0-Zr) 


400 [-6] 


406 (21) 


385 


-81 -9] 


-72 (20) 


-92 


-29 -3] 


-32 (22) 


-55 


PZT (R3m) 




















Ol(Zr-O-Ti) 


394 [-7] 


401 (0) 


401 


-90 [-6] 


-84 (23) 


-107 


-3[1] 


-2(1) 


-3 


02(Ti-0-Zr) 


466 [-1] 


467 (11) 


456 


-127 [-13] 


-114 (21) 


-135 


-2[1] 


-1 (-1) 


-1 



"This cluster-GGA calculation used a smaller Pb cc-pVDZ basis, 
rather than cc-pVTZ, due to convergence difficulty with the larger 
basis. 



exchange-correlation. Results from PBC-GIPAW with 
GGA exchange-correlation are also shown. Chemical 
shifts were determined using the corresponding values of 
aJJ^f in Table Iml The GGA 21 QM atom embedded 
cluster and PBC-GIPAW RMS errors differ by 4 ppm. 
As seen in Tables Hill the 21 QM atom embedded cluster 



chemical shifts are also in excellent agreement with those 
from the larger 65 QM atom cluster. 

PMN PBC-GIPAW calculations were based on a low 
symmetry (3x2x2) 60-atom perovskite supercell struc- 
ture with relaxed internal coordinates,— with B-site 
cations arranged according to the random-site modeli^ 



X-ray patterns of wcll-anncalcd PMN samples indicate a 
homogeneous average structure, which exhibits rocksah- 
hke 1:1 B-site ordering, which is well described by the 
random-site model. Locally, the random-site model cor- 
responds to B-site [111] planes, alternating between pure 
Nb and mixed Nb/Mg layers. Thus there are twice as 
many Nb-O-Mg as Nb-O-Nb 0-sites in the 60-atom su- 
percell, since the Mg atoms all reside in the mixed [111] 
planes. Similarly, there are no Mg-O-Mg 0-sites. Un- 
like earlier models, the random-site model satisfies charge 
neutrality locally. The previously accepted space-charge 
model was based on the apparent inability to fully anneal 
samplesi^ The 60-atom perovskite supercell structure)^ 
used in the present calculations, is consistent with the 
random-site model. Grinberg et ai— found good agree- 
ment between this relaxed 60-atom supercell structure 
and pair distribution functions (PDFs) obtained by neu- 
tron scattering experiments.— This indicates that this 
structural model reasonably represents the local struc- 
ture in PMN. Calculations for PMN were done only with 
PBC-GIPAW, and the results are shown in Table |Vl 



IV. DISCUSSION 

As previously noted in Ref. 0, there is a large 
anisotropy between the <7x,y and CTz principal values in 
Tabic ini The ax.y principal values are large and nega- 
tive (deshiclded), while (j^ values tend to be considerably 
smaller and positive (shielded). (As mentioned, the "z" 
principal axis is identified as that most closely aligned 
with the B-O-B bond direction of the target O atom.) 
The present calculations show that this anisotropy is also 
found in PZT and PMN and is reflected in the large Sax 
values (Eq.Hl) in Tables |IV] and El 

As shown in Ref. [Sj for prototypical perovskites, the 
large Sax values are due to paramagnetic contributions 
to ax.y from virtual transitions between 0(2p) and unoc- 
cupied B(?7d) states. The p-d hybridization contributes 
predominantly to (Tx.y, due to O atoms having only two 
nearest neighbors in perovskites, with linearly arranged 
B-O-B structural units. As the B-O-B bond distances 
vary, large variations in the chemical shielding can oc- 
cur, resulting in a strong dependence on rg. We find a 
similar dependence for the PZT and PMN alloy systems. 



A. PZT 

Figure [T] plots PZT isotropic and axial chemical shifts 
as a function of Tg, the shortest B-0 bond length of 
the targeted O atom. The B3LYP 21 QM atom clus- 
ter results from Table IIVI are plotted. For comparison, 
results for ST, BT, PT, and PZ and from experiment 
are also shown.— The dashed straight lines are the lin- 
ear fits to the calculated results for these prototypical 
perovskites, taken from Ref. i. A plot of PBC-GIPAW 
results from Table ITVl (not shown) is very similar, which 



is consistent with the generally good agreement of the 
chemical shift results in the Table between the two cal- 
culational approaches. A nearly linear dependence on 
Tg (with slope ~ 850 ppni/A) is seen in both Si^o and 
(5ax, across all the systems studied. The axial shift is 
plotted as 2(5ax ~ Sz — Siso (Eq. [2]) to emphasize that the 
linear dependence is largely due to Sx.y, while Sz has a 
much weaker dependence on rg, as previously noted for 
the prototypical perovskites^ As seen in the figure, the 
calculated PZT 50/50 results follow the same trends as 
in Ref. i. 

Pb(Zri_3;Ti3;)03 NMR ^^O magic angle spinning 
(MAS) central peak spectra were presented by Bald- 
win et alA for a range of concentrations x. Tabu- 
lated chemical shifts were given only for the endpoint 
PT and PZ compounds (reproduced here in Table HIT)) , 
whose spectra consist of well-defined narrow peaks. MAS 
removes broadening due to chemical shift anisotropy 
in powder samples, but only partly averages second- 
order quadrupolar broadening. The narrow peaks in 
the endpoint PT and PZ compounds indicate that elec- 
tric field gradients (EFGs) at the ^^O nuclei are small<^ 
This is consistent with first-principles calculations of O 
EFGsii^ Indeed, the "'^''O peak positions are within a 
few ppm of the experimental isotropic chemical shiftsi^ 
There are two inequivalent O atoms in PT, two "equa- 
torial" Ooq (which has two equidistant nn Ti atom) and 
one "axial" Oax atom (which has one short and one long 
Ti-0 bond). A narrow peak at 443 ppm has twice the 
(integrated) relative intensity of the peak at 647 ppm, 
and these were assigned to the Ocq and Oax atoms, re- 
spectively. These are in good agreement with the cal- 
culated results in Table IIIIl The B3LYP calculation ac- 
curately reproduces this splitting, while GGA underesti- 
mates it by ~ 45 ppm. In PZ, there are five inequivalent 
O sites, which corresponds to five MAS peaks, centered 
near 350 ppm, within ± ^ 40 ppm. Measured and calcu- 
lated values in Table Hill are in good agreement. 

At intermediate Ti and Zr concentrations, the spectra 
in Fig. 3 of Ref. Q show that only a few of the narrow 
endpoint peaks persist. As Zr is added to PT, the narrow 
647 ppm PT peak decreases quickly in intensity. It can 
no longer be observed in the x = 0.55 sample. A broad 
feature, between about 350 and 450 ppm is fully devel- 
oped near x ~ 0.50, with narrower embedded features at 
~ 370 and 430 ppm. This broad feature distribution of 
inequivalent 0-sites in the disordered PZT solid solution 
samples. The 430 ppm feature, which is close to the PT 
Ooq 447 ppm peak, is observed to persist down to 25% Ti 
concentration. Baldwin et alA assign the 430 ppm fea- 
ture to a site similar to that of the PT Oeq atom, i.e. a 
locally Ti-Ooq-Ti (undimerized) chain-like configuration. 
They further suggest that this peak could indicate Ti 
clustering on a spatial scale of at least two unit cells in 
PZT. 

Our results suggest an alternative explanation for the 
persistence of the observed 430 ppm feature. Our calcu- 
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FIG. 1: (Color online) Calculated PZT oxygen isotropic (5iso 
and axial 2(5ax chemical shifts (B3LYP 21 QM atom cluster 
values from Tables lrV)l , plotted as a a function of r^ , the short- 
est B-O bond length of the targeted O atom. For comparison, 
calculated results for ST, BT, PT, and PZ and from experi- 
ment (hollow symbols) are also shown.— The dashed straight 
lines are linear fits to calculated values for ST, BT, PT, and 
PZ from Ref H. 



FIG. 2: (Color online) GIPAW calculated PMN oxygen 
isotropic Siao and axial 2(5ax chemical shifts, plotted as a func- 
tion of rs- The PMN chemical shifts are from Table IVl and 
the symbols indicate the B-O-B' configuration, using the con- 
vention in that table. For comparison, results for cubic and 
tetragonal KNbOs are also shown. The dashed lines are the 
same as that in Fig. [l] As plotted, the rg values have been in- 
creased by 0.12 A (see text). The unshifted values are shown 
in the inset. 



lations show similar chemical shifts for a range of B-O-B' 
environments, with rg ranging between '^ 2.0 — 2.1 A in 
all the PZT 50/50 structural models, as seen in Table ITVl 
and Fig. [TJ Thus Ti clustering need not be invoked to 
explain the persistence of the 430 ppm feature in the 
measured spectra. 

The apparent disappearance, at intermediate concen- 
trations, of the 647 ppm peak does, however, indicate a 
reduced occurrence of a PT-like Oax site with a short 
1.77 A Ti-0 bond. We find such a site only in the relaxed 
P4mm and P2mm PZT 50/50 simulations. The P4mm 
and P2mm models have [100] type B-site ordering. The 
absence of the 647 ppm peak in the measurements indi- 
cates that local occurrences of [100] type B-site ordering 
are rare. Instead, the lack of the 647 ppm peak in our 
R3m PZT 50/50 simulations suggests that local rocksalt- 
like B-site ordering is more prevalent in disordered PZT. 
This conclusion is also supported by the R3m structural 
model having the lowest total energy of all the struc- 
tural models by ^ 23 mRyJ- Moreover, the calculated Ti 
EFG's for 50/50 R3m were significantly smaller and in 
better agreement with measured values, than the other 
B-site orderingsji 

Baldwin et alA also remark a narrow 287 ppm peak 
that appears in their x — 0.48 sample. This peak is 
broader in a; = 0.55 and x = 0.25 samples and is absent 
in PZ. These authors note that x = 0.48 is the compo- 
sition corresponding to the morphotropic phase bound- 
ary (MPB) and conjecture that the 287 ppm peak evi- 
dences a new oxygen environment in a distinct crystalline 
monoclinic phase, which has been suggested to bridge 
the MPBji They suggest that the new environment at 



X ~ 0.48 is a Ti-O-Zr site, which becomes ordered in 
the crystallographic sense as the long-range order of the 
monoclinic phase is established. As indicated by Fig. [1] 
isotropic chemicals shifts near 287 ppm are associated 
with large rg ~ 2.15 A. As most Zr-0 bond lengths are 
larger than those of Ti-0 in the PZ and PZT 50/50 
structural models, the 287 ppm feature could also be at- 
tributed to Zr-O-Zr sites with bond lengths distributed 
near ~ 2.15 A. 



B. PMN 

PMN, by contrast with PZT, is a heterovalent 2:1 mix- 
ture of the Nb^"'' transition metal cation and the Mg^+ 
alkaline earth cation. Figure [2] plots Pb(Mgi/3Nb2/3)03 
(PMN) isotropic and axial chemical shifts (GIPAW re- 
sults from Table W\ as a function of rg, the shortest B-0 
bond length of the targeted O atom. As plotted, the 
PMN Ts values, have been increased by 0.12 A for both 
5iso and (5ax, as further discussed below. The unshifted 
values arc plotted in the inset. The dashed lines are the 
same as that in Fig. [1] With the 0.12 A rigid shift for 
all 0-sites, the PMN chemical shifts are seen to follow 
the same linear trend as in the homovalent B-site sys- 
tems in Fig. [TJ For Nb-O-Mg coordinated 0-atoms, Tg is 
taken as the Nb-0 bond length, since the hybridization 
mechanism does not apply to Mg, which has no low-lying 
unoccupied d-states. (As mentioned, the linear depen- 
dence is due to paramagnetic contributions to (5xy from 
virtual transitions between 0(2p) and unoccupied B(77d) 



states.) Wc arc not aware of published NMR ^''O spectra 
for PMN. 

Although the linear dependence of Siso and (5ax in Fig. [2] 
have the same slope as the homovalent B-site systems, a 
0.12 A Tg offset is needed for the PMN values to fall on the 
same line. As mentioned, the linear dependence reflects 
variations in the magnitude of the paramagnetic 0(2p)- 
Nb(4d) hybridization contributions to the Sx,y principal 
values. This indicates that the effective Tg is controlled 
by the spatial extent of the paramagnetic screening cur- 
rents. The Nb^"*" cation could be expected to modify the 
spatial extent of these currents, compared to the B''+ 
cations in the homovalent systems, due to the larger elec- 
trostatic attraction of the Nb^+ cation. The 0.12 A off- 
set in PMN renormalizes, in effect, the strength of the 
0(2p)-Nb(4d) hybridization. This observation would ap- 
pear to indicate that for these systems the rg-dependence 
of dx^y is given by Sx,y = rn{rs — rp), where the slope m 
is nearly the same for all transition-metal coordinated 
O atoms in perovskites, while the intercept rg depends 
on other factors, such as the ionic charge of the nearest 
neighbor cation. To further examine this, we performed 
GIPAW calculations for cubic and tetragonal KNbOsr^ 
which are also plotted in Fig. [2l While the tetragonal 
axial O values (short Nb-0 rg) are consistent with the 
linear trend, the KNbOa tetragonal equatorial and cu- 
bic O results (two equidistant Nb-0 bonds) show signif- 
icant deviations. [Similar deviations occur in Fig. [1] for 
the tetragonal equatorial and cubic O BT results. We 
note that both K and Ba have larger crystalline ionic 
Shannon^ radii, 1.78 and 1.75 A, respectively, than Sr 
or Pb, 1.58 and 1.63 A, respectively.] Smaller deviations 
are also seen in the two largest rg values for the PMN 
oxygen atoms with two nearly equidistant Nb atoms, but 
for somewhat larger values of Vg than in KNbOs. In the 
ease of heterovalent B-site perovskites, the covalency of 
A-site atoms could also be important, as indicated by 
the smaller deviations, for long Vg, in PMN compared to 
KNbOs. These observations warrant further investiga- 
tion to clarify these issues. 

PMN is an end-point of the solid-solution series 
(x)PbTi03-(l-a;)Pb(Mgi/3Nb2/3)03 (PMN-PT). In 
the 60-atom PMN structural model, there are no Nb- 
O bonds as small as the short Ti-Oax bond ~ 1.7A in 
PT, which is associated with its high degree of tetrag- 
onality {c/a = 1.065). This short bond corresponds to 
the large 6iso — 640 ppm, which is also seen in the P4mm 
and P2mm PZT structural models, both of which also 
show a high degree of tetragonality. The tetragonality 
of PMN-PT decreases as the Ti concentration is reduced 
from PT-rieh compositions, and the average symmetry 
switches from tetragonal to rhombohedral at the mor- 



photropie phase boundary (MPB) x ~ 0.35. The largest 
piezoelectric response is typically achieved at concentra- 
tions near the MPB 4^ Polarization rotation has been 
proposed as the origin of the large piezoelectric response 
at the MPB, via intermediate monoelinic phases;^"— 
where the tetragonality increases as the polarization ro- 
tates from [111] to [100] (pseudocubic) directions, with 
applied electric field along a pseudocubic axis. Increased 
tetragonality, compared to PT, has been seen in some 
other perovskite based solid solutions, such as some Bi 
based materials . ^^'^^ The present calculations indicate 
that ^"^O NMR chemical shift measurements could be 
a useful probe in this regard, as increased tetragonal- 
ity is accompanied by shortened transition-metal/oxygen 
bonds. 



V. SUMMARY 

First-principles oxygen NMR chemical shift tensors 
were calculated for PZT and PMN, which are repre- 
sentative, respectively, of homovalent and heterovalent 
perovskite-structure B-site alloys. Quantum chemistry 
methods for embedded clusters and the GIPAW method 
for periodic boundary conditions were used. Results 
from both methods are in good agreement for PZT and 
prototypical perovskites. PMN results were obtained 
using only GIPAW. Both isotropic and axial chemical 
shift components were found to vary approximately lin- 
early as a function of the nearest-distance transition- 
metal/oxygen bond length, Tg. Using these results, we ar- 
gue against Ti clustering in PZT, as conjectured from re- 
cent -'^^O NMR measurements. Our findings indicate that 
^^O NMR measurements, coupled with first-principles 
calculations, can be an important probe of local struc- 
ture in complex perovskite solid solutions. 
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TABLE V: GIPAW calculated oxygen isotropic, axial and 
anisotropic components (ppm) of the chemical shift tensor for 
PMN. A 60-atom supercell with relaxed internal coordinates 
was usedi^ The B-O— B' bond notation is the same as in 
Table IIVI bond lengths are categorized as equidistant if they 
differ by less than 0.05 A. 



Nb-0 


B'-O 


(5iso 


(5ax 


^aniso 


A 


A 








Nb-O-Mg 










2.11 


2.09 


235 


-6 


-5 


2.11 


2.09 


235 


-7 


-2 


2.06 


2.08 


280 


-26 


-47 


2.07 


2.09 


277 


-22 


-43 


Nb-O-Mg 










1.89 


2.06 


419 


-109 


-27 


1.90 


2.07 


405 


-97 


-11 


1.90 


2.09 


390 


-110 


-25 


1.90 


2.11 


387 


-110 


-31 


1.91 


2.11 


390 


-111 


-83 


1.93 


2.08 


362 


-91 


-7 


1.93 


2.09 


368 


-94 


-10 


1.93 


2.10 


374 


-100 


-83 


1.94 


2.13 


369 


-75 


-36 


1.94 


2.18 


390 


-75 


-17 


1.95 


2.11 


362 


-81 


-11 


1.95 


2.18 


385 


-67 


-18 


1.96 


2.05 


346 


-78 


-29 


1.96 


2.12 


348 


-61 


-44 


1.97 


2.05 


338 


-71 


-36 


1.97 


2.10 


343 


-64 


-11 


1.97 


2.13 


327 


-53 


-13 


1.98 


2.06 


340 


-51 


-47 


1.99 


2.13 


316 


-43 


-7 


2.01 


2.06 


310 


-36 


-39 


Nb-O-Nb 










1.86 


2.23 


436 


-130 


-21 


1.87 


2.14 


421 


-130 


-3 


1.87 


2.14 


421 


-129 


-8 


1.87 


2.23 


437 


-130 


-18 


1.89 


2.20 


416 


-114 


-27 


1.90 


2.20 


414 


-110 


-24 


1.91 


2.11 


394 


-111 


-23 


1.91 


2.11 


392 


-108 


-18 


1.94 


2.03 


401 


-122 


-8 


1.95 


2.01 


399 


-122 


-11 


Nb-O-Nb 










2.06 


2.09 


347 


-62 


-9 


2.08 


2.08 


340 


-56 


-11 



